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Abstract. A simple model of ballistic aggregation and fragmentation is proposed. 
The model is characterized by two energy thresholds, E^gg and -Bfrag, which demarcate 
different types of impacts: If the kinetic energy of the relative motion of a colliding 
pair is smaller than i?agg or larger than ii^frag, particles respectively merge or break; 
otherwise they rebound. We assume that particles are formed from monomers 
which cannot split any further and that in a collision-induced fragmentation the 
larger particle splits into two fragments. We start from the Boltzmann equation 
for the mass- velocity distribution function and derive Smoluchowski-like equations for 
concentrations of particles of different mass. We analyze these equations analytically, 
solve them numerically and perform Monte Carlo simulations. When aggregation and 
fragmentation energy thresholds do not depend on the masses of the colliding particles, 
the model becomes analytically tractable. In this case we show the emergence of the two 
types of behavior: the regime of unlimited cluster growth arises when fragmentation is 
(relatively) weak and the relaxation towards a steady state occurs when fragmentation 
prevails. In a model with mass-dependent i?agg and iJfrag the evolution with a cross- 
over from one of the regimes to another has been detected. 
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1. Introduction 

Collision-induced aggregation and fragmentation are ubiquitous processes underlying 
numerous natural phenomena. For a gentle collision with a small relative velocity, 
colliding particles can merge; a violent collision with a large relative velocity can cause 
fragmentation. For intermediate relative velocities, particles usually rebound. These 
collisions may still be irreversible — the kinetic energy could be lost in inelastic collisions. 
Important examples of such systems are dust agglomerates in the Earth atmosphere 
or in interstellar dust clouds and proto-planetary discs [1, 2, 3, 4]. Another example 
is dynamic ephemeral bodies in planetary rings, see e.g. [5, 6, 7]. A comprehensive 
description of the aggregation and fragmentation kinetics in such systems is very 
complicated. Therefore it is desirable to develop idealized models that involves three 
types of coUisions in the simplest possible way. 

The understanding of the ballistic-controlled reactions is still quite incomplete [8] . 
Ballistic aggregation has attracted most attention (see [3, 4, 9, 10, 11, 12, 13, 14] 
and references therein) and a few studies were also devoted to ballistic fragmentation 
(see [15, 16, 17, 18]). The situation with aggregation and fragmentation operating 
simultaneously has been analyzed only in a very special case when all particles have 
the same relative velocity and the after-collision fragment mass distribution obeys a 
power-law [7]. Moreover, studies of pure ballistic fragmentation arc usually based on the 
assumption that all particles may split independently on their mass and relative velocity 
between the colliding grains [16, 17]. In reality, the type of an impact strongly depends 
on the relative velocity [4]; furthermore, the agglomerates are comprised of primary 
particles ("grains") that cannot spht into smaller fragments [5, 6]. The fragmentation 
model with the splitting probability depending on energy has been studied in [19]; this 
model, however, does not consider ballistic impacts of many particles, but rather an 
abstract process of a successive fragmentation of one body with a random distribution 
of the bulk energy between fragments. 

In this paper we propose a model of ballistic aggregation and fragmentation which 
accounts for three types of collisions depending on masses and relative velocity of a 
colliding pair. In section 2 we introduce the model, write the Boltzmann equation for 
the joint mass-velocity distribution function, and deduce from the Boltzmann equation 
the rate equations for concentrations of various mass species. Section 3 is devoted to 
the theoretical analysis of rate equations. Numerical verification of theoretical results 
and simulation results in situations intractable theoretically is given in section 4. The 
last section 5 concludes the paper. 

2. The model 

Consider a system comprised of primary particles (monomers) of mass mi and radius ri, 
which aggregate to form clusters of 2, 3, . . . , A;, . . . monomers with masses = kmi. In 
some applications (e.g. in modeling of dynamic ephemeral bodies [5, 6]) it is appropriate 
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to consider clusters as objects with fractal dimension D; for compact clusters D = 3. 
The characteristic radius of an agglomerate containing k monomers scales with mass 
as Tfc ~ rik^/^ . We assume that when the kinetic energy of two colliding clusters in 
the center-of-mass reference frame (the "relative kinetic energy" in short) is less than 
E'agg, they merge. In this case a particle of mass {i + j)mi is formed. If the relative 
kinetic energy is larger than ^'agg; but smaller than E'^ag, the colliding particles rebound 
without any change of their properties. Finally, if the relative kinetic energy exceeds 
-Sfrag, one of the particles (we assume that the larger one) splits into two fragments. We 
denote Pi^k-i the probability that a particle of mass k splits into particles of masses i 
and k — i. Obviously, J2iPi,k-i = 1 and Pi^k-i = if /c < i. 

We restrict ourselves to dilute and spatially uniform systems. Let fi = f{vi,t) 
be the mass-velocity distribution function which gives the concentration of particles of 
mass rui with the velocity Vi at time t. The mass-velocity distribution function evolves 
according to the Boltzmann equation 

where /^^^, 11^^^ and f^^^ are respectively the coUision integrals describing collisions 
leading to aggregation, rebound, and fragmentation. The first integral reads 

i+j=k 

X fi (vi) fj (vj) © (-Bagg - Eij) S{mkVk - rUiVi - mjVj) 

~ 5Z ^fcj / ^^i / ^^0 {-Vkj ■e)\vkj-e \ x 

j 

X fk (4) {Vj) © (£^agg - Ekj) . 

Here aij — ri (i^^^ + j^^^^ is the sum of radii of the two clusters, while ruk — mi + mj 
and mkVk = VifUi + mjVj, due to the conservation of mass and momentum. Wc also 
introduce the relative velocity, Vij — Vi — Vj, the reduced mass, /x^j = miUij/irni + m^), 
and the relative kinetic energy, E^j = ^fiijvfj. The unit vector e specifies the direction 
of the inter-center vector at the collision instant. The factors in the integrand in Eq. (2) 
have their usual meaning (see e.g. [20]): afj \vij ■ e \ defines the volume of the collision 
cylinder, © (— % ■ e ) selects only approaching particles and © (^'agg ~ Eij) guarantees 
that the relative kinetic energy does not exceed i^agg to cause the aggregation. The 
first sum in the right-hand side of Eq. (2) refers to collisions where a cluster of mass k 
is formed from smaller clusters of masses i and j, while the second sum describes the 
collisions of /c-clusters with all other aggregates. 
For collisions leading to fragmentation wc have 

^k^^i^k) = Y.Y.Pk,j-k (l - crl j dvj Jdvi JdeO {-Vij ■ e) x 

j i<j 

X 1% • e I fj (vj) fi (vi) © {Eij - £;frag) ^{vi, Vj, Vk) (3) 



i<k 
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X fk {Vk) fi (Vi) © {Eki - -Bfrag) , 

where ruj — + ruj^k s-nd we use the abbreviation, A{vi, Vj, Vk) — 5{mjVj + rUiVi — 
mkVk — nT'k-jVk_j + mivl) for the factor which guarantees the momentum conservation at 
the colUsion. The after-colhsional velocities tT^_^- and t7/ are determined by a particular 
fragmentation model. The first sum in Eq. (3) describes the collision of particles of mass 
i and j (j > k, j > i) with the relative kinetic energy above the fragmentation threshold 
-E-frag- The larger particle, i.e. the particle of mass j, splits with the probability Pkj~k 
into two particles of mass k and j — k thereby giving rise to a particle of mass k. The 
second sum describes the opposite process, when particles of mass k break in collisions 
with smaller particles. In the present study we do not need an explicit expression for 
the velocities vl,_j and -tT/ of the fragments. We do not also need an expression for the 
collision integral I^^^; it has the usual form (see e.g. [20]) with a slight modification to 
account for the requirement that the relative kinetic energy Eij belongs to the interval 

(-^agg ^ Eij < -£'frag)- 

Thus we have a mixture of particles of different masses and each species generally 
has its own temperature. For this (granular) mixture we write 

rii ^ dvifi{vi) , N^^rii (4) 

i 

where rii is the number density (concentration) of particles of mass i and N is the total 
number density. Using the mean kinetic energy of different species one can also define 
the partial granular temperatures Tj for clusters of mass i and effective temperature T 
of the mixture [21]. We assume that the distribution function fi{vi,t) may be written 
as [4, 14, 21] 

M^i^ ^) = 3 , Ci = —, (5) 

where Vq .i{t) = 2Ti{t)/mi is the thermal velocity and (j){ci) the reduced distribution 
function. For the force-free granular mixtures the velocity distribution functions of the 
components are not far from the Maxwellian distribution [21], which reads in terms of 
the reduced velocity c = ?7/t'r, 

0M(c)=7r-3/2exp(-c2). (6) 

The equipartition between different components may, however, break down, in the sense 
that the partial temperatures Tj are not equal and differ from the effective temperature 
T [21]. Here we ignore the deviation from the Maxwellian distribution \ and possible 
violation of the equipartition and use the approximation, 0j(cj) ^ 0m (q) and Ti ^ T 
for all i. We also assume that the temperature of the system does not depend on 
time. This is formally inconsistent within the Boltzmann equation (1), yet in many 

I Note that for a slightly niodified model, where only in a small fraction of collisions that fulfil 
the aggregation criterion, particles merge and only in a small fraction of collisions that fulfil the 
fragmentation criterion, particles split, the velocity distribution is close to the Maxwellian, like in 
a granular mixture. At the same time this modified model would lead to the same kinetic equations 
(9), but with the renormalized time scale. 
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applications the temperature is approximately constant on a astronomical time scale, 
e.g. this happens in planetary rings where the viscous heating due to the shearing mode 
of the particle orbital motion keeps the granular temperature constant [22, 4]. The 
consistent approach would be to modify the Boltzmann equation to take into account 
gradients of the local hydrodynamic velocity, which will result in additional terms in 
the velocity distributions /j, proportional to these gradients. If we assume that such 
gradients are very small but still sufficient to support constant temperature due to 
viscous heating, we can neglect the small corrections to the distribution functions and 
approximate them with a gradient- free form (5). 

Integrating Eq. (1) over Vk we obtain the equations for the zero-order moments 
of the velocity distribution functions f^, that is, for the concentrations n^. Taking 
into account that collisions resulting in rebounds do not change the concentrations of 
different species and using (2)~(3) together with (4)-(6) we arrive at rate equations 

d 1 °° ooi /1\ 

'll'^k = 9 H CijUiUj -UkY^ Ck,ini + XI XI ^ij^i^j ( 1 - ^^ij ) Pk,j-k 



i+j=k i=l j=k+li=l 

1 

— ( 

2 



- nfc(l - Sik) J2 ^k,i^i ( 1 - \^i,k ] , (7) 



with rates given by 




It is useful to verify that the above kinetic equation (7) fulfills the condition of mass 
conservation, J2k kmiUk — M — const., where M is the total mass density. 

The probability of splitting pik depends on geometric and mechanical properties 
of the aggregates and generally it is quite complicated. For concreteness we focus on 
splitting into (almost) equal fragments. Namely, we assume that a particle of mass 2kmi 
splits into two equal halves, while a particle of mass {2k + l)mi splits into particles of 
mass kmi and (A;-|-l)mi. For this choice of the splitting probability, the kinetic equation 
reads 

d 1. 00 k 

^ i+j=k 1=1 i=l 

2k 2k+l 

+ 2^A-2k,in2kni (1 - (^2fe,i/2) + ^ ^2k+i,in2k+ini (1 - Wi,i/2) 
1=1 1=1 

2fe-l 

+ X Mk-i,in2k-ini (1 - (^2fe-i,i/2) . (9) 

i=l 

In the next sections we study this model theoretically and numerically. 
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3. Theoretical analysis 

To understand the qualitative behavior it is instructive to start with the simplest model 
which allows an analytical treatment. 

3.1. Constant rates 

Consider first the model with constant rates Aij and Cjj. Without loss of generality 
we can choose these rates to be 

Cij = 2 , Ai^j = 2A , (10) 

The parameter A quantifies the relative intensity of fragmentation with respect to 
aggregation. Fragmentation prevails when A > 1 while aggregation wins in the opposite 
case of A < 1. If A = 1 two processes are in a balance. 

Even in this simple case we still ought to analyze a cumbersome system of infinitely 
many equations. To gain insight it is useful to consider the evolution of the total density 
N — X^njfc. (In many problems involving aggregation and fragmentation this quantity 
satisfies a simple equation that does not contain other densities.) Summing up all 
equations (9) we obtain 

AT = -(1 - A)Ar2 - An? 

which has indeed a neat form, although it additionally involves the density of monomers. 
This density evolves according to 



hi = -2niN + 2A 



712 (2ni + ri2) + ns{ni + 712) + ^nj 



(12) 



Although we do not have a closed system we can already reach some qualitative 
conclusions. Equation (11) indicates that two regimes are possible. If A < 1, i.e. 
when aggregation prevails, the system continues to evolve leading to formation of larger 
and larger clusters; when A > 1, one expects that the system reaches a steady state. We 
now analyze these situations in more detail. 

3.1.1. Unlimited cluster growth, A < 1 In this case larger and larger clusters will arise. 
Since the total mass is conserved, one expects that the concentration of small clusters 
will rapidly decrease. Therefore, rii <^ N when i ^ 1 and therefore one can omit the 
second term on the right-hand side of (11). Similarly one can keep only the first term 
on the right-hand side of (12). This leads to the simphfied equations 

TV ~ -(1 - A)A^^ hi ~ -2niN (13) 

which are solved to yield the large time behavior: 

n, ~ r"/<'-^' (15) 
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Further, one anticipates that the density distribution approaches the scahng form 
nfc~r'^$(a;), x^^ (16) 

in the scahng hmit t oo, k — > oo, with the scaled mass x = k/t^ kept finite. (Here z 
is the dynamic exponent characterizing the average mass: (k) ~ t^.) The scahng form 
agrees with mass conservation: Y^kn^ — J dxx^{x) is manifestly time-independent. 

The exponent z can be found from the known asymptotic behavior of N{t). Indeed, 
writing 

POO 

AT = V rife ~ dx $(a;) - 

and matching this with already known asymptotic behavior (14) we conclude that z — 1. 
If we further assume that $(x) ~ x'' for x <^ 1 and combine this asymptotic with 
z — 1 and the scahng ansatz (16) we obtain rii ~ t~'^t~^. Matching with (15) we get 
// = 2A/(1- A). Therefore 

.,^2A/(l-A) 

- ? (l) (^^) 

when k <^ t. Obviously, the above equation implies the asymptotic time dependence 
rik ~ i^VCi--^) and the mass dependence rik ~ ior x — k/t <^ 1. 

3.1.2. Relaxation to a steady state, A > 1 For A > 1 the system evolves to a steady 
state with constant concentration of clusters. In this case hk — N = and Eq. (11) 
yields, 

m = NVl - A-i (18) 

The densities rik rapidly decay with k. Therefore n2k <^ rik ior k ^ 1 and the governing 
equations (9) for the stationary concentrations simplify to 

fc-i 

J2 riiUk-i - 2(1 + X)nkN = (19) 

i=l 

where we have approximated a finite sum up to A; ^ 1 by an infinite sum and ignore the 
terms containing n2fe,n2fe±i- The above equation is supposed to be valid for large k; it 
is certainly invalid for k = 1 when the right-hand side does not vanish. The qualitative 
form of the large k asymptotic behavior is determined by the mathematical structure 
of (19). To extract this asymptotic let us consider the simplest version when Eq. (19) 
is valid for all /c > 2. Specifically, let us probe the model 

k-l 

J2 riiUk-i - 2(1 + X)nkN = -(1 + 2X)N5k,i , (20) 

i=l 

where the amplitude (1 + 2A) was chosen to set A'^ = 1. (For model (20), this choice 
merely sets the overall amplitude.) 
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The infinite system (20) forms a recurrence and therefore it is solvable. Introducing 
the generating function 



X(z) = 5^ rife, 



k>l 



we recast (20) into a quadratic equation 

Af^ - 2(1 + X)Af + (1 + 2X)z = 
which is solved to yield 

AT = (1 + A) - ^{l + Xy-{l + 2X)z 
Expanding ^{z) we arrive at 



1 + A 



47r 



1 - 



1 



;i + A)^ 



(21) 



(22) 



(23) 



(24) 



T{k + 1) 

From this solution one gets ni = (A + l/2)/(A + 1), which of course directly follows from 
Eq. (20) as well. For large k, equation (24) simplifies to 

1 k 



rik ~ ^^^^ k 



47r 



1 



We considered other tractable versions when Eq. (19) is exact above a certain 
threshold, k > k, + 1, while for k = 1, . . . , k we use the same modification as in Eq. (20) 
for /c = 1. In this case instead of (22) one gets A/"^ - 2(1 + X)Af + P{z) = with 
P{z) = AiZ + . . . + Af^z'^. The root of P{z) = (1 + A)^ closest to the origin is positive [one 
can show that it exceeds unity, z^ > 1] and non-degenerate. Expanding the generating 
function J\f = 1 + A — ^^{1 + Xy — P{z) leads to the asymptotic ~ k~^/'^z~'^. The 
above argument favors the asymptotic behavior 

Ak-^'l^e-^^ 



(26) 



This asymptotic form is universal and only the parameters A, 7 depend on the specificity 
of the model, that is on the parameter A. 

It is impossible to determine A, 7 since models like (20) are uncontrolled 
approximations. Let us still use such models and choose the simplest one which obeys 
the exact relation of Eq. (18). The model (20) is inappropriate as it fails to satisfy (18): 
(A + l/2)/(A + 1) > VI - A-i. Modifying Eq. (19) at A; = 1, 2 yields (we still set = 1) 

(27) 



with 



^ riiUj - 2(1 + X)nk = -q5k,i - (2A + 1 - q)5k,2 

i+j=k 



g = 2(l + A)\/l- A-1 



ensuring the validity of (18). The same approach as before gives (26) with A and 7. In 
particular, 



7 = In 



(1 + A) 



yi-A-i + A- yi-A-^ 

A 



(28) 
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where we have used the short-hand notation A = 1 + 2A — 2(1 + A)-\/l — For A = 2 
(which we have studied numerically) one gets 

7 = 0.495156... (29) 

This is an uncontrolled approximation, of course. Interestingly, the result is rather close 
to the numerically obtained value 7 ^ 0.465. 

3.2. Mass-independent energy thresholds 

We now turn to the analysis of the situation when aggregation and fragmentation energy 
thresholds E'agg and E'^ag are constant. In this case the total density of clusters evolves 
according to 

-1 00 CXD 1 

^ = -0(1 - >^)EEC^,Jn^n, - -\C,,,nl , (30) 
^ i=i j=i ^ 

where A~^ = ^^^b-Jt (^1 _ (^1 _^ Es_^jT) e~^^sg/^). This again imphes the existence of the 
two opposite evolution regimes: For A > 1 the relaxation to a steady state is expected, 
while for A < 1 - the regime of the unlimited cluster growth. 

The rates Cjj = C{i,j) and Aij — A{i,j) differ by a constant factor A; moreover, 
they are homogeneous functions of their arguments: 

A(ai, aj) = a'^Aii, 3) , C(ai, aj) = a'^C(i, 3) (31) 

with the exponent 

-^-i. (32) 

which follows from the relation for a particles mass rui — imi, the cross-section of the 
collision cylinder, af^ ~ (i^/^ + i^/^^ , and Eqs. (8). Plugging the scahng Ansatz (16) 
into Eq. (9), taking into account that for /c ^ 1 the summation may be approximated 
by integration, and exploiting the homogeneity of the rate kernels, Eq. (31), we obtain 
(see e.g. [23, 25, 24] for analysis of similar integro-differential equations for the re-scaled 
mass distribution) 

^ (2$(x) + x^'ix)) = dy^y) [(C.,, + A.,,)$(x) (33) 

-\Cy,x-y^{x -y)- AA2:c,v^{2x) 



(34) 



From Eq. (33) we find the scaling exponent 
1 _ /3 _ 2\-^ 

^ ~ i-v ~ V2 ~ 

Although the scaling theory does not allow to determine the scaling function $(a;), one 
can find the total concentration of clusters from Eq. (34): 

N{t) ~ t-' ~ ^-2D/(3i?-4) _ (35) 

Correspondingly, the average clusters mass grows as {k) — M/N ~ t^. 
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3.3. Dependence of the energy thresholds on masses of colliding particles 

In the preceding analysis we have assumed that i^agg and -Efrag do not depend on the 
mass of coUiding particles so that Xij — Aij/Cij = A is constant. In reality, however, 
such dependence does exist, implying that X^j is a function of i and j. Still, if Ajj > 1 or 
Xij < 1 for all i and j, the qualitative behavior of a system is similar to that for the case 
of constant i^agg and -Efrag^ For X^j > 1 a relaxation to a steady state is expected, while 
for Xij < 1 an unlimited cluster growth is observed. The most interesting behavior is 
expected when Xij — 1 changes its sign with increasing clusters masses i and j. In this 
case one anticipates a cross-over from one type of evolution to another. 

To choose a realistic dependence of E'agg and E'frag on the masses of colliding 
particles, one needs more details of the coUision process. We shall use the threshold 
energy for baUistic aggregation that takes into account surface adhesion [26]. In this 
case 



where E^^^ is expressed in terms of the monomer radius, particle surface tension, the 
Young modulus and the Poisson ratio of the particle material (see [26] for the exphcit 
expression for E^gg)- 

For the energy of fragmentation we assume that it is equal to the energy required 
to create an additional surface, which may be roughly estimated as twice the area of 
the equatorial cross-section of the larger particle (recall that the model assumes, that 
the larger particle in a collision pair breaks down). Hence we adopt the following mass 
dependence for -Efrag^ 



where % = 1 if i > j, = if i < j and On — 1/2; E^^^^ — 2'kJs1^i, with 7^ being the 
surface tension. 

4. Numerical simulations 

In our numerical studies we apply two different approaches - the solutions of the 
system of differential equations and the direct modeling of random aggregation and 
fragmentation processes (with the corresponding rates Cjj and Aij) by means of Monte 
Carlo (MC) method. In the former case we use 1000 equations and in the later one 
100, 000 monomers (we always used the mono-disperse initial conditions). The approach 
based on the solution of differential equation has an obvious deficiency as one must 
approximate an infinite system of equations with a finite one. The MC approach is more 
time consuming, yet it has an advantage of directly imitating the physical processes in 
which particles are involved. To model the fragmentation and aggregation kinetics by 
MC we use the standard Gillespie algorithm [27, 28] (see [29] for the apphcation of this 
algorithm to the aggregation and fragmentation processes). 



4 




(36) 




(37) 
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The results presented on Figs. 1-11 confirm our theoretical predictions qualitatively 
and quantitatively. For constant rates two opposite types of evolution have indeed been 
observed: the relaxation to a steady state for dominating fragmentation (A > 1) and the 
unlimited cluster growth when aggregation prevails (A < 1). Both numerical approaches 
(the solution of the differential equations and MC) yield very close results. 

In Fig. 1 the evolution of the concentration of clusters of different mass is shown 
for the A < 1 regime when the cluster growth continues ad infinitum. Note that all 
concentrations nk{t), except for ni{t) which always decays, initially increase and then 
decay to zero. Figure 2 shows that the decay of cluster density N{t) well agrees with 




Figure 1. Evolution of cluster concentrations rik in the case of constant rates with 
A = 0.5. The lines correspond to the numerical solution of 1000 differential equations, 

symbols arc the results of MC simulation with 100 000 monomers for the monodispcrse 
initial conditions. Note, that while ni{t) always decays, nfe(t) initially increase and 
than decay to zero. 

the theoretical prediction (14). Figures 3, 4 show respectively the asymptotic evolution 
of cluster concentrations and the distribution of the cluster mass for x — k/t <^ 1. 
The theoretical predictions, Eqs. (15) and (17), are in a good agreement with the 
simulations. Relaxation to a steady state in the case when fragmentation dominates 
(A > 1) is illustrated in Fig. 5, while Fig. 6 demonstrates the corresponding stationary 
cluster mass distribution. 

Note that the numerical simulations confirm the theoretical form of the steady state 
cluster mass distribution. Qualitatively similar behavior is observed for the case of the 
baUistic kinetic coefficients, Eqs. (8), with the constant aggregation and fragmentation 
energies. Again, for A < 1, as for constant kinetic coefficients, clusters unlimitedly grow. 
Fig. 7, while for A > 1 the system relaxes to a steady state. Fig. 8. The cluster mass 
distribution in the steady state may be anew very well fitted with the nearly-exponential 
form, Eq. (26), see Fig. 9. 

In Fig. 10 the prediction (35) of the scaling theory is compared with the numerical 
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t 



Figure 2. Evolution of the total number of clusters N for the same system as in 
Fig. 1. The solid line corresponds to the numerical solution of the differential equations, 
symbols - to the MC simulation and the dotted line shows the theoretical prediction, 
Eq. (14), N{t) ~ for t > 1. 




iE-6-^ — . . — . \ I 

100 1000 

t 



Figure 3. The long-time limit behavior of cluster concentrations for the same 
system as in Fig. 1. In accordance with the theoretical predictions, Eqs. (15) and (17), 
the cluster concentrations nk{t) decay for t 1 with the same slope t~2/(i-A) ^ ^-4^ 
shown by the solid line. 
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lE-S-^ 

1 10 100 1000 

k 

Figure 4. The cluster mass distribution at different time instants for the same system 
as in Fig. 1. The initial cluster distribution at t = 10 (long-dashed line) drastically 
differs from that in the scaling regime, t ^ 1. The dottcd-dashcd and dashed lines 
show respectively the cluster mass distribution for t — 500 and t — 1000. The solid 
line shows the theoretical prediction, Eq. (17), rife k^, with /x = 2A/(l — A) = 2. 
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10 



100 



Figure 5. Evolution of cluster concentrations nk{t) for the case of constant kinetic 
coefficients with A = 2. After a certain period of time the system relaxes to a 
steady state. Lines correspond to the numerical solution of 1000 differential equations, 
symbols - to the results of MC simulation (100 000 monomers) for the monodisperse 
initial conditions. 
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Figure 6. The steady state distribution of cluster mass for the same system as in 
Fig. 5. The sohd line corresponds to the numerical solution, the dotted line - to the 

theoretical prediction, Eq. (26), Uk = Ak~^/'^e~'^^ . The constant 7 — 0.465, obtained 
by fitting, is very close to the theoretical value of 7 = 0.495, Eq. (29). 

data. Again we see that the agreement between the theory and simulations is rather 
satisfactory. Finally Fig. 11 illustrates evolution of the system with the ballistic 




Figure 7. Evolution of cluster concentrations nk(t) for the case of ballistic kinetic 
coefficients with constant aggregation and fragmentation energies E^^^jT = 0.9, 
Eivag/T = 3 and A < 1. The cluster dimension is D = 3. 

coefficients that depend on the cluster mass in accordance with Eqs. (36) and (37). It is 
interesting to note that the system tends initially to a quasi-steady state, as previously 
for the case of A > 1, but then a cross-over to a different evolution regime, corresponding 
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to A < 1 takes place. In the latter regime all cluster concentrations decay with a similar 
slope, close to t~^, still to be explained theoretically. 
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Figure 8. Evolution of cluster concentrations nk{t) in the case of ballistic kinetic 
coefficients with constant aggregation and fragmentation energies Eg^g^/T = 0.3, 
Eiiag/T = 3 and A > 1. The cluster dimension is f = 3. Similar to the case of 
constant kinetic coefficients, the system relaxes to a steady state. Lines correspond to 
the numerical solution of f 000 equations, symbols - to the results of MC simulation 
(100 000 monomers) for the monodisperse initial conditions. 




Figure 9. steady state cluster mass distribution of for the same system as in Fig. 
8. Similar to the case of constant kinetic coefficients with A > 1, the distribution nk 
has a nearly-exponential form, Eq. (26), nj, = Ak~^/^e~'''' . The constant 7 = 0.315 is 
obtained by fitting. 
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Figure 10. Evolution of total number of clusters N(t) for the case of ballistic 
kinetic coefHcients with constant aggregation and fragmentation energies E^^^/T = 0.9, 
Eha,g/T = 3 and A < 1 for different cluster dimensions D. Solid lines correspond to 
the numerical solution of 1000 differential equations, dotted lines show the prediction 
of the scaling theory, Eq. (35). 
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Figure 11. Evolution of the cluster concentrations nk{t) for the case of ballistic 
coefficients C^, Aij with the mass-dependent aggregation and fragmentation energies 
E,^{i,j) and i;frag(i,.7), given by Eqs. (36) and (37) with E^JT = 0.1, E^JT = 
0.6. Lines correspond to the numerical solution of 1000 differential equations, symbols 
- to the results of MC simulation (100 000 monomers) for the monodisperse initial 
conditions. The cluster dimension is = 3. Note that the system tends initially to a 
steady state, clearly seen for the first few clusters masses. Later its evolution alters to 
the regime corresponding to the unlimited cluster growth. In this regime the cluster 
concentrations decay with a slope close to ~ , showed in the figure by the dotted 
line). 
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5. Conclusion 

We analyzed the dynamics of a system where particles move ballistically and undergo 
collisions which can lead to decrease or increase of the number of particles. The precise 
outcome depends on the kinetic energy E'kin in the center-of-mass reference frame. We 
proposed a simple model with two threshold energies, E'agg and E'frag, which define a type 
of an impact: For E'kin < E^agg the colliding particles merge, for E'agg < £^kin < -E'frag 
they rebound, and for -Efrag < -E'kin one the particles (the larger one) splits upon the 
collision. We assume that the aggregates are composed of 1, 2, . . . , A;, . . . monomers and 
split into two equal (for an even number of monomers in the cluster) or almost equal 
(for an odd number of monomers) pieces. The monomers are assumed to be stable, that 
is, they do not further split. For this model we wrote the Boltzmann kinetic equation 
for the mass- velocity distribution function of the aggregates and derived rate equations 
for the time evolution of the cluster concentrations. The ballistic rates were obtained 
in terms of the aggregation and fragmentation energy thresholds £^agg and £^frag, masses 
of the colliding particles and the temperature of the system (which was assumed to be 
constant). The MaxweUian velocity distribution for all species in the system was also 
assumed. 

Wc analyzed theoretically and studied numerically the rate equations. In the 
numerical studies wc used two different methods - the solution of the system of 
differential equations and Monte Carlo modeling. Both numerical methods yielded 
very close results. We started with the simplest case of constant rates and observed 
two opposite evolution regimes — the regime of unlimited cluster growth and of the 
relaxation to a steady state; we described both these cases analytically. For the regime 
of the unlimited cluster growth we obtained the asymptotic time dependence for the 
cluster concentrations and for their mass distribution. For the relaxation regime, which 
corresponds to the prevailing fragmentation, we derived the asymptotic behavior of the 
stationary mass distribution. In the evolving regime, the cluster concentrations decay 
as a power law in time; the stationary mass distribution has a nearly exponential form. 
Theoretical predictions are in a good agreement with numerical results. 

We also studied the case of mass-dependent rates arising in the situation when 
aggregation and fragmentation energy thresholds are constant. We observed that the 
behavior of the system is qualitatively similar to that of the system with the constant 
rates. Surprisingly, we detected that the steady state cluster mass distribution has also 
a near-exponential form. We developed a scaling theory for the asymptotic large-time 
behavior of the cluster concentrations and checked it numerically for different fractal 
dimensions of the aggregates. The numerical data agree well with the results of our 
theory. 

Finally, we explored numerically the case of the baUistic kinetic coefficients with the 
aggregation and fragmentation energies depending on the mass of colliding particles. For 
the aggregation energy threshold we use the available in literature result for a collision 
of particles with surface adhesion. For the fragmentation energy threshold we adopted a 
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model where -Efrag is proportional to the surface energy of the maximal cross-section of 
the larger particle in the colliding pair. For this model the dependence on mass of £^frag 
is much stronger than that of -Eagg- As the result, the evolution of the system, where 
the fragmentation initially prevails and drives it to a steady state, alters at later time 
when the unlimited cluster growth eventually wins and then it continues ad infinitum. 
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